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Kinetic energy of solid neon by Monte Carlo 
with improved Trotter- and finite-size extrapolation 
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The kinetic energy of solid neon is calculated by a path-integral Monte Carlo approach with a 
refined Trotter- and finite-size extrapolation. These accurate data present significant quantum 
effects up to temperature T = 20 K. They confirm previous simulations and are consistent with 
recent experiments. 



In a previous workO we reported theoretical results 
about the average kinetic energy of rare gas solids (kryp- 
ton, argon and neon), modeled by a Lennard- Jones 
(LJ) interaction. For heavier crystals the thermodynam- 
ics was .-approached by means of the effective-potential 
methodcla. This approach allows us the use of all clas- 
sical methods through the construction of an approxi- 
mate effective classical phase-space distribution (see for 
details^). Monte Carlo simulationsM'Effil with the ef- 
fective potential were favourably compared with path- 
integral Monte. Carlo (PIMC) simulations, and also ap- 
plied to argonB, reproducing very well the experimental 
density and specific heata. 

In spite of their large mass, krypton and argon show 
relevant quantum effects at easily accessible tempera- 
tures; for instance, the average kinetic energy is much 
larger than its corresponding classical value. Our calcula- 
tions of the average kinetic energy of argon suggested the 
realization of neutron Compton scattering (NCS) exper- 
iments, whose outcomes have been later found in perfect 
agreement with our predictionstj. For increasing value 
of the quantum coupling, anharmonic second order cor- 
rections arising from the odd part of the potential were 
found relevanfJJ and have been recently inserted in the 
effective potential formalismM. 

For neon, where thepStrong anharmonicity shows up 
even in the ground statetj, we preferred to resort to path- 
integral Monte Carlo IPIMC) simulations, with a refined 
Trotter extrapolation!!!; the procedure consists in adding 
to the PIMC data the contribution of the harmonic ap- 
proximation after the subtraction of the corresponding 
finite Trotter data. 

The available experimental data for the kinetic energy 
of solid neon obtained by NCS experimentsEj, were in 
evident disagreement with our PIMC results; all theoret- 
ical data were lower than the experimental ones. The 
reason_.was not clear at all. Since the experiments with 
argonEj were not yet performed, we supposed that the 
LJ potential were an over-simplification of the true po- 



tential and that many-body interactions could play an 
important role. ._. 

In a recent work Timms et al$3 report new NCS mea- 
surements of the kinetic enewy of solid neon, which dif- 
fer from the previous onestia Both experiments were 
carried out and analyzed in the regime of the so-called 
impulse approximation, which becomes exact when the 
energy and the momentum transferred to the sample 
are infinite. Timms et al. reached higher momentum 
and energy transfers, and they claim that this improve- 
ment has made the final-state-effect corrections to the 
observed longitudinal Compton profile irrelevant: hence, 
the data analysis had less sources of uncertainty and was 
thus more reliable. They also did new PIMC simula- 
tions, using different potentials (Aziz and LJ) in order to 
definitively establish whether the original disagreement 
between theory and experiment was due to a too rough 
model potential: it turned out that the computed ki- 
netic energy depends very weakly on the model poten- 
tial. Their experimental and PIMC data are consistent 
between themselves, and also confirm the validity of our 
previous simulations— 

In a recent paperlij we have developed a systematic 
method for improving the Trotter number extrapolation 
of PIMC data; in addition, we have recently extended this 
procedure in order to take into account also the effect of 
the finite size of the simulation boxO, so that we are now 
able to obtain much more accurate results. 

The PIMC method is based upon the semi-group prop- 
erty of the density matrix 

p{q',q;(3) = J dq P _ 1 ...dq 1 p(q' ,q P _ x ;T)...p{q 1 ,q;r), 

(1) 

where /3 = 1/T and r = [3/P, P being the Trotter num- 
ber. In order to make the above formula of practical use, 
the density matrix element p(q e , <jf^_i; t) is usually taken 
in the lowest high-temperature approximation, giving 
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rise to the so-called primitive action. The latter tends to 
the exact density matrix as P — > oo. This is the formal- 
ism we are dealing with here. In order to get values of the 
averages of physical observables, many simulations at dif- 
ferent values of P must be performed, and then the data 
must be extrapolated as P — ► oo. The more "quantum" 
is the system, the larger must be the maximum value 
of P. The finite P estimates G(P) of the averages can 
be expanded as G(P) = G(oo) + g 2 /P 2 + gi/P 4 + ■ • -O, 
and frequently the term in 1/P 2 is not sufficient for a 
satisfactory fit. ■— ■ 

The method we suggestedEJ overcomes this problem. 
The idea is to take advantage of the fact that the ther- 
modynamics of a harmonic system can be obtained ana- 
lytically, even at finite P (see Ref. (L8|); nevertheless, such 
a system shows a very strong dependence on P, and the 
results obtained at finite P will not be close to those at 
P = oo unless the condition P 3> / = phuj/2 is fulfilled 
for any system's frequency w. However, at low tempera- 
tures, when quantum effects are most important (in the 
temperature region where one should use the highest val- 
ues of P), the harmonic approximation (HA) of a solid 
system is surely meaningful, though rough. Indeed, as 
T — > 0, harmonic excitations play a very important role 
in the thermodynamics. As long as the "quantum char- 
acter" of the system increases, this beconaes-less and less 
true, but the self-consistent HA (SCHA)lii>H2l eventually 
allows to recover a simple harmonic-like system whose 
behaviour is very similar to that of the real system. 

Our idea is to improve the extrapolation as P — > oo 
in PIMC simulations accounting for the P dependence 
of the harmonic contributions to the PIMC estimates of 
physical observables. The procedure consists in adding 
to the rough PIMC data G{P) the deviation from the 
P = oo estimates calculated for the SCHA of the system: 



G SC (P) = G(P)+ GfM°°)~Gf6(.P) 



(2) 



In such a way, the improved estimates Gsc(P) will show 
a much weaker dependence on the Trotter number P, the 
scaling behaviour in 1 / P 2 is reached earlier and the maxi- 
mum Trotter number necessary to get the correct asymp- 
totic result is lower. We remember that as P increases 
the computer time grows both because of the larger sim- 
ulation box and of the worse statistics. 

Another important point which has not yet been 
deeply investigated in relation to quantum simulations 
is the dependence of the data on the simulation box size. 
It is well known that for systems undergoing a phase 
transition the finite size of the simulated sample has dra- 
matic effects, because in the critical region the correlation 
length diverges, and in order to simulate such a system 
particular procedures known as "finite-size scaling" must 
be used. In the classical case, if the system is far from a 
phase transition, the problem of how to reach the ther- 
modynamic limit regime is in general easily overcome, 
without using enormous samples (size effects, with mod- 
ern computer capability are in general not a problem). 



However, dealing with quantum systems, subtle phenom- 
ena can occur: for instance, in d-dimensional lattices the 
discreteness of the Brillouin zone due to the finite par- 
ticle number introduces a nonphysical gap into the dis- 
persion curve which gives rise, at low temperature, to an 
exponential behaviour of the specific heat, instead of the 
correct T d scaling Bloch law, obtained for a linear dis- 
persion of the soft modes. This effect can be observed in 
simulations!^. ._. 

Following the idea of Eq. (|J) we suggests to correct 
the raw PIMC data, at the SCHA level, also with re- 
spect to their dependence on N, the number of particles 
composing the actually simulated sample, 



G SC (P, N) = G(P, N) + (oo, oo) - Gf>(P, N) 



(3) 

Preliminary testsEl made on I-d nonlinear systems con- 
firm that even with a chain composed by very few parti- 
cles it is possible to get the thermodynamic limit of the 
averages of observables only adding to the raw simulation 
data this harmonic-like correction, practically making the 
extrapolation as N — > oo unnecessary. 

To model solid neon, we considered a fee lattice com- 
posed by N particles (labeled with 3-d indices i, j) inter- 
acting through a pairwise potential, 



V r (9) = 5E w (l«i-*D 



(4) 



Several choices, as we know, ace possible for the model 
potential. Since Timms et alx3 showed that the depen- 
dence of the kinetic energy on the potential is not critical, 
we chose to model our system by the LJ 12-6 potential 
v(r) = 4e[((i/r) 12 — (c/r) 6 ] with the potential parame- 
ters e and a taken from the literatureEl (e = 36.68 K 
and a = 2.787 A). We neglect the dynamic effect of 
the interactions beyond nearest-neighbours, whose con- 
tribution to the potential energy is taken into account by 
a static-lattice approximation. We performed constant- 
density simulations evaluating the pressure within each 
run. The density was adjusted in such a way to get a 
practically vanishing pressure, (the pressure is always less 
than 0.07p* ss 15 atm, being p* — e/a 3 the characteris- 
tic pressure) in order to best reproduce the experimental 
settings; the zero-pressure densities turned out to be very 
close to the experimental ones. The sample was an fee 
lattice of 108 atoms with periodic boundary conditions; 
in order to test the above described finite-size correction 
scheme we made test runs changing the box size up to 864 
particles. We used the Metropolis algorithm, with both 
single- and many-particle moves. The maximum Trotter 
number P was 48. Each run consisted of 200, 000 steps 
per particle for equilibration followed by 1,200,000 fur- 
ther steps, during which the averages were accumulated 
every 5 steps. For each run we estimated the statisti- 
cal uncertainty, taking into account the MC correlation 
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times; these vary with P and N and never exceed 400 
steps. 

In order to make finite- Trotter and finite-size harmonic 
corrections in the spirit described above, we need the 
SCHA potential, 
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where £ ; Q = qf — q^ ; and qo,i is the equilibrium position 
of the i — th particle, which is fixed, being determined by 
the particle density, z = 12 is the coordination number 
and w and Q 2< if are adjustable parameters determined 
imposing that the average of the actual potential V(q), 
and of its first and second derivatives are equal to the 
corresponding averages obtained for V (q). all averages 
are performed using the finite- P density distribution cor- 
responding to Vo(q), which is a gaussian. A shorthand 
way of expressing these gaussian averages as a formal 
power series turns out to be very useful in this case, 
namely (/(£ ld )) = exp(£> Q/3 dadp/2) /(0) (summation 
over repeated indices) where £ ;d = £i +d ~ £i (d labels the 
nearest-neighbour displacements) and V a(i — (£j d £^) 
is the variance matrix of the gaussian distribution. 

Since v(r) — > oo as r approaches 0, averages like 
(v(|d+£ id |)) would diverge. This is an artifact of 
the harmonic approximation: the gaussian is small but 
nonzero at r = 0, where the true distribution would van- 
ish; the formal power series exp(V a/3 d a dp /2)V (q) is then 
only asymptotic. However, the nonphysical contributions 
from the potential core can be simply eliminated by trun- 
cating the series. Hence we can expand the averaged po- 
tential and its derivatives up to second order in the T>'s, 
finally obtaining the following SCHA equations 
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where v = v(d), v' = v'(d) and so on; d = |d| is the 
nearest-neighbour distance. The indices k (wavevector) 
and v (polarization) label the normal modes, obtained 

by diagonalizing Sl 2 ^ to the eigenfrequencies Q 2 ^- In 
particular, 4i/£ = $Z d [l — cos(k ■ d)] and A£ results from 
the polarization diagonalization. The quantityLJ't£l 



h coth(P/i£) 



2mfl^ cosh 



(8) 



is the normal coordinate mean square fluctuation at finite 
P, with sinh(/j£) = fiMV^j (2P) ; the corresponding limit 
for P — > oo is easily recovered. Moreover, 
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2?|| and T>± are the mean square fluctuations of the com- 
ponents of £ id , parallel and orthogonal to d, respectively. 
Due to the symmetry properties of the fee lattice, the 108 
2?'s reduce indeed to three only, and making a further 
isotropy approximation we assume that the two trans- 
verse components are equal, 2?j_,i — D12 — P_l- 
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The SCHA finite— P estimates can be obtained as loga- 
rithmic derivatives of the partition function 



Z {P,N) = e -/w»™/3 Y[ [ 2 sinh(P^)]- 



(12) 



The well-known partition function for a system of quan- 
tum harmonic oscillators is recovered as P — » 00. In 
order to get the finite-A values we use the discrete 
mesh in the Brillouin zone corresponding to that value 
of A; the thermodynamic limit is obtained by the spe- 
cial points method. The SCHA kinetic energy is K, = 

(2N)- 1 £ kj/ mfi 2 kQk where «k = a k( F > N )- In this wa y, 
we are able to get both finite P and finite N corrections. 
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FIG. 1. Reduced kinetic energy K/e vs. particle num- 
ber N. The Trotter number is P = 8, the temperature 
is T = 20 K and the reduced density is p/p* — 0.945. 
p* = rn/a 3 = 1.5479 gem -3 is the characteristic density. The 
triangles are raw PIMC data, the squares are PIMC data plus 
finite- Trotter corrections (Q) and the circles are PIMC data 
plus finite- Trotter and finite-size corrections (^). Error bars, 
when not shown, lie inside the symbols. Lines are guides for 
the eye. 

Through this analysis we have concluded that for the 
kinetic energy of neon a simulation box with 108 atoms 
and periodic boundary conditions is barely large enough 
to mimic the thermodynamic limit behaviour (i.e. the 
finite size corrections are of the order of the statistical 
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error), as it can be seen in Fig. [l] where it's shown how 
the finite-size corrections scheme works. The relative ef- 
fect of finite-size and finite- Trotter corrections is well seen 
also in Fig. [^, where it is shown how the finite- Trotter 
corrections work. The small difference between the ex- 
trapolations at P = oo represents the size effect. 

The results of our simulations are shown in Fig. |^: they 
are consistent with the experimental daiall and confirm 
the validity of other PIMC simulationsLxlif . 

In conclusion, the SCHA correction scheme gives an 
estimate on how large both the finite-size and the finite- 
Trotter effects are. It is indeed important to control how 
much the data are affected by the finitcncss of the simula- 
tion box, especially for observables which possibly show 
a stronger iV-dependence than the kinetic energy, such 
as the specific heat or the frequency moments, or for sys- 
tems with larger quantum coupling. 

We gratefully acknowledge useful discussions with P. 
Nielaba. 
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FIG. 2. Reduced kinetic energy K/e vs. the inverse square 
of the Trotter number P. The simulation box contains 108 
particles, the temperature is T — 5 K and the reduced den- 
sity is pj p* = 0.968. The triangles are raw PIMC data, and 
the circles are PIMC data plus finite Trotter and finite size 
corrections (|^). Error bars, when not shown, lie inside the 
symbols. The lines are fits of corresponding data. 
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FIG. 3. Kinetic energy K, of solid neo«-*s. temperature T. 
The crosses are experimental data. fromLJ, the open squares 
are PIMC simulations data fromcj and the triangles are our 
PIMC data with refined Trotter extrapolation. Error bars, 
when not shown, lie inside the symbols. The line is a guide 
for the eye. 
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